%% segement all the files

FileList = {
   
    'Ext_Annulaire.mat'
    'Ext_Auriculaire.mat'
    'Ext_IMRL.mat'
    'Ext_Index.mat'
    'Ext_Majeur.mat'
    'Ext_Poignet.mat'
    'Ext_Pouce.mat'
    'Fermeture_Poing.mat'
    
    % prolematic
    'Flex_Annulaire.mat'
    'Flex_Auriculaire.mat'
    'Flex_IMRL.mat'
    'Flex_Index.mat'
    'Flex_Majeur.mat'
    % end of problematic
    
    'Flex_Poignet.mat'
    'Flex_Pouce.mat'
    'Incl_Radiale.mat'
    'Incl_Ulnaire.mat'
    'Ouverture_Poing.mat'
    'Pronation.mat'
    'Supination.mat'
    
     %'Fichier_Test.mat'
     %'Cocontraction.mat'
    };

Fs = 2000;

window = Fs;
X = []; % contains the trials from all classes
Y = []; % contains the classification from all classes

FileList = {'Flex_Pouce.mat' 'Flex_IMRL.mat' 'Flex_IMRL.mat'};

for curFileIx = 1:length(FileList)
    %load(['D:/dropbox/Gipsa-work/data_EMG/franck/' FileList{curFileIx}]);
    load(['D:/dropbox/Gipsa-work/data_EMG/marc/' FileList{curFileIx}]);
    
    datar = data(:,1:12);
    
    [dataf] = filter_emg(datar,Fs);
    
    if (size(data,2) == 13)
        trigger_channel_available = true;
        trigger_channel = data(:,13:13);
    else
        trigger_channel_available = false;
    end;
    
    if (trigger_channel_available)
        [best_index,best_trigger] = generate_trials_postitions_using_continous_trigger_channel(trigger_channel);
    else
        [best_index,best_trigger] = generate_trials_postitions_auto(dataf);
    end;

    plot_trials(data,FileList{curFileIx},best_trigger);

    % cut signal
    %index = find(diff(trigger)==1); % because we do it earlier now!!!!
    trials_per_class = zeros(size(dataf,2),window,length(best_index));
    disp(['Trials: ' int2str(length(best_index)) ' ' FileList{curFileIx}])
 
    %get the data for these triggers
    for i=1:length(best_index)
        %calculate a range around the trigger
        range = 100;
        trials_per_class(:,:,i) = dataf(best_index(i)-range:best_index(i)-range+window-1,:)';
    end
    
    X = cat(3,X,trials_per_class); % put all trials together
    Y = cat(1,Y,ones(length(best_index),1) * curFileIx);
end

%% Classification over X and Y

% covariance estimation
COV = covariances(X);


%shuffuling
NTrial = size(X,3); % trials equals movements
ix = randperm(NTrial);

% permutate the data for classification
COV = COV(:,:,ix); % the sequence is updated from 1 ... ntrial to a random permutation
Ys = Y(ix); % it takes Y in the permutated form

% cross validation
NCV = 10;
% MDM
disp('---------- MDM --------------');
out = cross_valid(COV,Ys,NCV,@mdm); % COV is produced from X and Ys is produced from Y
disp('Confusion matrix :');
disp(confusionmat(Ys,out));
disp(['Accuracy : ' num2str(100*mean(Ys==out)) ' %']);


% FGMDM
disp('---------- FGMDM --------------');
out = cross_valid(COV,Ys,NCV,@fgmdm);
disp('Confusion matrix :');
disp(confusionmat(Ys,out));
disp(['Accuracy : ' num2str(100*mean(Ys==out)) ' %']);

% K-NN
disp('---------- KNN --------------');
% 10 neigbourg
out = cross_valid(COV,Ys,NCV,@knn,10);
disp('Confusion matrix :');
disp(confusionmat(Ys,out));
disp(['Accuracy : ' num2str(100*mean(Ys==out)) ' %']);

% LDA on std
disp('---------- LDA on std (no Riemannian geometry) --------------');
Feat = sqrt(squeeze(sum(X(:,:,ix).^2,2)))';
%  Constants
NTrial = size(X,3);
NTrialTest = fix(NTrial/NCV);
out = zeros(size(X,3),1);
%cross validation
for iCV = 1:NCV 

    idxTest = 1+(iCV-1)*NTrialTest:iCV*NTrialTest;
    if iCV ==NCV
        idxTest = 1+(iCV-1)*NTrialTest:NTrial;
    end
    idxTraining = 1:NTrial;
    idxTraining(idxTest)=[];

    out(idxTest) = classify(Feat(idxTest,:),Feat(idxTraining,:),Ys(idxTraining));

end

disp('Confusion matrix :');
disp(confusionmat(Ys,out));
disp(['Accuracy : ' num2str(100*mean(Ys==out)) ' %']);

%% MDM Accuracy Vs. window size
return;

Wsize = [0.05:0.05:0.8]*Fs;
offset = 0;
NTrial = size(X,3);
ix = randperm(NTrial);
for i=1:length(Wsize)
    
    % covariance estimation
    COV = covariances(X(:,fix(offset+(1:Wsize(i))),:));


    %shuffuling

    COV = COV(:,:,ix);
    Ys = Y(ix);

    % cross validation
    NCV = 10;
    % MDM
    out = cross_valid(COV,Ys,NCV,@mdm);
    acc(i) = (100*mean(Ys==out));
    
end
figure;
plot(1000*Wsize./Fs,acc);
xlabel('Window size (ms)');
ylabel('Accuracy');
box on;grid on;


%% LDA Accuracy Vs. window size


Wsize = [0.05:0.05:0.8]*Fs;
offset = 0;
NTrial = size(X,3);
ix = randperm(NTrial);
for i=1:length(Wsize)
    
    Ys = Y(ix);

    % cross validation
    NCV = 10;
    % MDM

    % LDA on std
    Feat = sqrt(squeeze(sum(X(:,fix(offset+(1:Wsize(i))),ix).^2,2)))';
    %  Constants
    NTrial = size(X,3);
    NTrialTest = fix(NTrial/NCV);
    out = zeros(size(X,3),1);
    %cross validation
    for iCV = 1:NCV 

        idxTest = 1+(iCV-1)*NTrialTest:iCV*NTrialTest;
        if iCV ==NCV
            idxTest = 1+(iCV-1)*NTrialTest:NTrial;
        end
        idxTraining = 1:NTrial;
        idxTraining(idxTest)=[];

        out(idxTest) = classify(Feat(idxTest,:),Feat(idxTraining,:),Ys(idxTraining));

    end
    acc(i) = (100*mean(Ys==out));
end
figure;
plot(1000*Wsize./Fs,acc);
xlabel('Window size (ms)');
ylabel('Accuracy');
box on;grid on;